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Abstract 

We propose a novel stochastic metliod to generate patlis conditioned to start in an initial state and 
end in a given final state during a certain time tj. These paths are weighted with a probability given 
by the over damped Langevin dynamics. We show that these paths can be exactly generated by a 
non-local stochastic differential equation. In the limit of short times, we show that this complicated 
non-solvable equation can be simplified into an approximate stochastic differential equation. For 
longer times, the paths generated by this approximate equation can be reweighted to generate the 
correct statistics. In all cases, the paths generated by this equation are statistically independent 
and provide a representative sample of transition paths. In case the reaction takes place in a solvent 
(e.g. protein folding in water), the explicit solvent can be treated. The method is illustrated on 
the one-dimensional quartic oscillator. 
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I. INTRODUCTION 



The problem of finding the pathway of chemical or biological reactions is of utmost 
importance for the understanding of their underlying mechanisms, as it allows to have better 
control on these reactions [Ij. For instance, in the realm of proteins, understanding the 
pathway between the unfolded state and the native state, or between two native states of 
the protein (allostery) may help prevent certain reactions or on the contrary favor them. 
Recent progress in single molecule experiments have allowed to monitor the spontaneous 
thermal folding and unfolding of single proteins, or the force induced unfolding of proteins 
[M]. 

In the following, we will study the spontaneous or the driven transition between an initial 
state denoted A and a final state denoted B. 

This problem has been addressed mainly by stochastic methods which start from an 
initial path and deform it by sampling the vicinity of the path. These are the path sampling 
methods |3H7]. The main drawback of these methods is that they are time consuming, 
and they generate strongly correlated trajectories. As a consequence, the space of sampled 
trajectories depends strongly on the initial used path. The same kind of problem exists for 
the Dominant Pathway method [HI |9], where the minimal action path depends strongly on 
the initial guess. 

From now on, we assume that the system is driven by stochastic dynamics in the form of 
an overdamped Langevin equation 

dx 1 dU , , 

*=-^&+"<" « 

For the sake of simplicity, we illustrate the method on a one-dimensional system, the gener- 
alization to higher dimensions or larger number of degrees of freedom being straightforward. 
In this equation, x{t) is the position of a point at time t in a potential U{x), 7 is the friction 
coefficient, related to the diffusion constant D through the relation D = ksT/'-f, where ks 
is the Boltzmann constant and T the temperature of the thermostat. In addition, ri{t) is a 
Gaussian white noise with moments given by 

< >= (2) 
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< vitHt') >= -^6{t - t') (3) 

7 

It is well known that the probability distribution P{x,t) for the particle to be at point x 
at time t is given by a Fokker-Planck equation [lOj 

dKT^p'^p) (4) 



dt dx \ dx dx 
where /3 = l/ksT is the inverse temperature. In this one dimensional model, the initial state 

A is characterized by its position xq at time and the final state B by its position xj aX time 

tf. This equation is thus to be supplemented by a boundary condition P(x, 0) = 5{x — Xq) 

where Xq is the initial position of the particle. 

It is convenient to go to the Schrodinger representation, by defining 

The function \E'(x,t) satisfies the imaginary time Schrodinger equation 
with 



Using the standard notations of quantum mechanics, one can conveniently write 

P{xf, tf\xo, 0) = e-^(^(-/)-^(-"))/2 < x/|e-*^^|xo > (7) 
where the Hamiltonian H is given by 

In eq.([7]), we have denoted by P(x f,tf\xQ,0) the probability for a particle to start at Xq at 
time and end at x/ at time tf, to emphasize the boundary conditions. 

It is well-known that the ground state of H, which has energy, is '^o{x) = e~^^^^'^^'^/ y/Z 
where Z is the partition function of the system, and all eigenstates \E'q, of H have strictly 
positive energies > 0. The spectral expansion of P can be written as 

g-/3C/(x) 

P(a;/,t/|a;o,0) = + ^ e-*^^'^P„(x/, Xq) 
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We see that for large tf the system converges to the Boltzmann distribution, and that its 
relaxation time is given by the inverse of the first eigenvalue = 1/Ei. In systems with 
high energy barriers, such as proteins, the gap Ei may be very small, and consequently the 
time Tn which in this case is identified with the folding time, can be very long. 

Using the Feynman path integral representation, we may thus write eq.Q as [TT] 

P{xf,tf\xo,0) = e-/5(f^(-/)-f^(-o))/2 f^'''''^ Vxit)exp ( --^ f'dt ( 7^2 + -^(x)^^ 

(9) 

In the following, we will be mostly interested in problems of energy or entropy barrier 
crossing, which are of utmost importance in many chemical, biochemical or biological reac- 
tions. As we already mentioned before, the archetype of such reactions is protein folding, 
a model we will use in the rest of this paper. A protein is a small biopolymer, which es- 
sentially may exist in two states, namely the native state (with biological activity) and the 
denatured state (with no biological activity) [12j. The protein being a small system (up 
to a few hundred amino-acids) , it never stays in one of the two states, but rather makes 
rare stochastic transitions between the two states (see figjl]). The picture which emerges is 
that of the system staying for a long time in one of the two states and then making a rapid 
transition to the other state. 



It follows that for most of the trajectory, the system makes uninteresting stochastic 
oscillations in the well, and can be described by normal mode analysis. Rarely, there is 
a very short but interesting physical phenomenon, which is the fast transition from one 
minimum to the other. 

This picture has been confirmed by single molecule experiments j2l |3], where the waiting 
time in one state can be measured, but the time for crossing from one state to the other is so 
short that it cannot be resolved. This scenario has also been confirmed recently by very long 
millisecond molecular dynamics simulations [13] which for the first time show spontaneous 
thermal folding-unfolding events. 

According to Kramers theory, the total transition time tk (waiting + crossing) scales like 
the exponential of the barrier energy 
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Figure 1: A long Langevin trajectory in the double- well. 

The "Kramers time" tk is the sum of two times: 

• the waiting time in the potential well 

• the crossing time over the barrier tq 

It is well known that the crossing time tq is small compared to Tk and indeed, Hummer |14| 
and subsequently Szabo [I5j have shown that 

Tc ~ In^-^ << Tk 

These Kramers and crossing times are averages. In fact, these times are distributed 
(random variables) and single molecule experiments or long molecular dynamics simulations 
allow to compute their probability distributions. 

However, it seems a bit wasteful to simulate proteins over huge time scales (milliseconds), 
during which only small conformational vibrations occur, just to observe interesting physical 
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crossing events which occur very rarely, on the sub-microsecond scale. 

The goal of this paper is to show how one can generate a representative sample of tran- 
sition paths, starting in state A at time and ending in state B at some arbitrary time tf. 
The typical times of interest are not the (long) folding times, but rather the (very short) 
transition or barrier crossing times. In mathematical terms, we are looking for the paths 
starting from A at time and conditioned to end in state B at time tf << tk- 

II. THE CONDITIONAL PROBABILITY 

Using the path integral representation of eq.(|9|, we see that the probability for a path 
{x{t)} starting at Xq at time 0, to end at at is given by 

Piixit)}) = le-/3(^(-/)-f^(-o))/2 exp ( ^ f' dt ( + -V(x)] ) (10) 

A \ AkeT Jo \ 7 J J 

where 

A= f (ia;/e-''(^("^)-^(""»/2 Vx{t) exp ( dt ( 7x2 + -^(x)^ ^ (11) 

J J{xo,o) \ 4:kBT Jq \ 7 J J 

The conditional probability over all paths starting at Xq at time and ending at xj at 
time tf, to find the system at point x at an intermediate time t is given by 

n^^t) = — — ^ — -Q(x,t)p(x,t) 

P(x/,t/|xo,0) 

where 

P(x,t) = P(x,t|xo,0) 

Q{x,t) = P{xf,tf\x,t) 
The equation satisfied by P is given by Q, whereas that for Q is given by 

ot ox^ ox ox 

It follows easily that the equation for the conditional probability V{x,t) is given by 
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dV d fdV d 

dt dx\dx^ dx ^ 

Comparing this equation with the initial Fokker-Planck Q and Langevin ([T]) equations, 
one sees that it can be obtained from a Langevin equation with a modified potential 

dx D dU ^d\nQ , ^ 
This equation has been previously obtained using the Doob transform |16j and is known 



in the probability literature as a Langevin bridge: the paths {x{t)} generated by (13) are 
conditioned to end at {xf,tf). It is the new term in the Langevin equation that guarantees 
that the trajectory starting at (a;o, 0) will end at {xf,tf). 
Using eq.Q for Q, one can write equation (13) as 



^ = 2^ A In < x.|e-(*/-*)^|x > +v{t) (14) 
dt J ox 

Using the analogous of the correspondence principle of quantum mechanics [T7], i.e. 
I — )■ p, this equation can also be rewritten in the form 

^ =< xit) > +vit) (15) 

where by definition 

1 / 1 /•*/ / . 1 

<x>= 1 ^ / Vx{T)x{t) exp -TTT^ / dr jx^ + -V{x) 



(16) 



Note that for large time tf, the matrix element in eq.(14) is dominated by the ground state 
of H, namely < x/|e~^*-f~*^^|x >~ e~^^^^^^'^^^^^^^ and as expected one recovers the standard 
(unconditioned) Langevin equation. 

Since we have a natural splitting of the Hamiltonian H as H = Hq + Vi with Hq = 
~^7^J&" ^^'-^ ^1 ~ y/'i'lksT, it is convenient to rewrite the above equation as 



2 — 7^ In < x,|e-(*/-*)^o|x > +2^- In . + r/(t) (17) 







X > 


< 




X > 



dt 7 dx 

Note that the first term in the r.h.s. above is singular at t = tf and is thus responsible 
for driving the system to {xf,tf) whereas the second one is regular. It follows that the first 
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term is the only term which can drive the system to {xf,tf), and any approximation which 
keeps the second term finite for t = tf will not affect this property. 

This nice bridge equation cannot be used "as is", since we don't know how to compute 
the function Q or equivalently the matrix element in the above equation. There are many 
ways to approximate this function. It is important however, to preserve detailed balance as 
well as possible, that the approximation retains the symmetry of the matrix element. 



III. THE MODIFIED LANGEVIN EQUATION AND REWEIGHTING 

The only approximation we found which remains local in time, i.e. which does not 
give rise to an integro-differential stochastic equation is the symmetric form of the Trotter 
approximation, commonly used in quantum mechanics [llj. Indeed, for short times t, a very 
simple and convenient symmetric approximation for Q is given by 

which translates into 



< Xf\e \x >~ e ** ^ ' s^y y i'^ \ >> 

It would be nice to relate the range of validity of this equation to the spectrum of H. 
Indeed, as was shown before, the spectrum of H corresponds to all the dynamical times of 
the system (folding times, transition times, etc.). We have not succeeded in finding such 
a relation except in the solvable case of the harmonic oscillator. In that case, it can easily 
be shawn that the natural expansion parameter is tA where A is the constant gap between 
the energy levels of H . As mentioned before, in the case of protein folding, the folding time 
which is the inverse of the first gap of the system can be very long, and we might expect the 
above approximation to be valid for times much smaller than this time. In particular, this 
approximation would allow to investigate the crossing times, much shorter than the folding 
time. 



Plugging eq.(18) in eq.(26) we obtain the approximate Langevin bridge equation which 
in arbitrary dimension (or with arbitrary number of degrees of freedom) reads 
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§-^-W'^'f-'m^m (19) 

where f](t) is a white noise vector whose components satisiy the relations Q and ([s]) and 

V{x) = {VUf - 2kBTV^U (20) 



The first term in the r.h.s of eq.(19) is the one which drives the particle to reach x/ at 
time tf. The potential which governs this bridge equation is not the original U{x) but rather 
the effective potential V{x). Note also that the force term is proportional to {tf — t) and 
thus becomes small as the particle gets close to its target site. 

In order to build a representative sample of paths starting at (xq, 0) and ending at {xf,tf), 
one must simply solve this equation for many different realizations of the random noise. Only 
the initial boundary condition is to be imposed, as the singular term in the equation imposes 
the correct final boundary condition. An important point to note is that all the trajectories 



generated by eq.(19) are statistically independent. From a numerical point of view, this 
means that this equation can be fully parallelized, and from a statistical point of view, it 
implies that all trajectories can be used in the representative sample. This last important 
point is to be contrasted with most existing methods where the sample are generated by 
some stochastic (Monte Carlo) methods which generate highly correlated trajectories. 

Before presenting examples of application of this method, let us discuss how to correct 
for the fact that the total time tf should be small for the approximation to be valid. 

Due to this restriction, the statistic of trajectories is not exact for larger times. Indeed, 



if eq.(19) were exact, each trajectory would be generated with its correct weight, and if 
one wanted to calculate observables, one would just have to compute simple white averages 
over all trajectories. However, as the equation is approximate, one needs to resample the 
ensemble of trajectories, that is, assign them a new weight. As we will show, the resampling 
weight is easily obtained. 



Indeed, if we consider the sample of trajectories generated using eq.(19) between (a;o,0) 



and (xfjtf), the weight of each trajectory should be given by eq.(lO). However, it is clear 



from eq.(19) that, using the Ito prescription, the weight with which it was generated is given 
by 



9 



exp 



7 



, f dx Xf — X tf — t 
dti ^ + ^ 



VV{x 




AksT Jo \ dt tf-t ' 
Up to a normalization, the reweighting factor for a trajectory is thus given by 



(21) 



exp 
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tf 



dt 



dx 
~dt 



Xf — X tf — t 



tf-t 



+ 



472 



VV{x 




(22) 



This quantity is easily calculated and allows for a correct evaluation of averages over 
paths. 

This reweighting technique can also be used to generate paths statistically exactly sam- 



pled according to eq.(13). Indeed, consider eq.(15). The expectation value < x{t) > can 
be computed by generating at each time t an ensemble of (approximate) trajectories start- 



ing from the current point x at time t and ending at xj at time t/ by using eq.(19). By 



reweighting them using the weights of eq.(22), we can reliably compute < x{t) > and thus 



solve eq. ( 15 ). Note that this procedure which generates correctly weighted trajectories might 



seem computationally costly. However, since all trajectories are independent, they can effi- 
ciently be generated using massive parallelization. 

IV. THE NATIVE STATE 



Eq.(18) is in fact not quite valid between non normalizable states like \x > and \xf >, in 
that it is not true to order 0{t^). However it is true between a normalizable state and |a; >. 
Assume that the final state of the system is defined by a probability distribution 0(x). For 
instance, for the case of a protein, (j){x) could represent the Boltzmann weight around the 
native state of the protein. The probability for the system to start at x at time t and end 
at time tf in the native state is given by 



Q{x,t) = / dy (l){y)P{y,tf\x,t) 



or usmg 







where restricts the integration over y to the vicinity of the native state. 



(23) 
(24) 



With this definition of Q, it is straightforward to see that eq.(12) and (13) are still valid 
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Using the approximation (18) we can write 



Q{x,t) = e/5^(-)/2e-i^i(-) f ^ 0(y)e-^^(^)/2e-i^^(^)e~'^^ (25) 



where A = ^/4:7^{tf - t)/(3'j. 



As the function (p restricts the integration in (25) to the vicinity of the native state, we 
can approximate the potential U{x) in this region by a quadratic expansion in terms of the 
normal modes 



U{x) ~ —{x — Xf) 



It follows that V and Vi are also quadratic and thus the integral (25) can be performed. 
Although we will consider only one-dimensional cases in the examples, we present the results 
for the mult i- dimensional case. 

Denoting by Uij = q^..^^. |^/ the Hessian matrix of normal modes around the native state, 
the potential U can be written in that region as 



U{x) = \ ^{Xi - x{)Uij{Xj - xj) 

and thus 
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V{x) = ^^(xj — x{)Qij{xj — Xj) — 2kBTTi Ui 



where the symbol Tr denotes the trace of the normal mode matrix and 



k 

The function Q can be easily calculated as 
where 



W., = J2D^kiI+^-^D),| 

k ^ 
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where / is the unit matrix and 



tf-t 
JJij — ujij H ilij 



The bridge equation becomes then 



rlr- 1 1 

^ = - E ^^^(4 - -4^2itf- t)V.V{x) + v^{t) (26) 



V. INCLUDING THE SOLVENT 

In many cases, in particular protein folding, one wants to include explicitly the solvent 
molecules, most often water. It is thus desirable to generate trajectories which are condi- 
tioned for the protein coordinates, but not for the water molecules. 

Denoting by Xj the set of coordinates of the water molecules, the conditional probability, 
over all paths starting at {xq, Xq} at time and ending at Xf at time tf (irrespective of the 
position of the solvent molecules), for the system to be at {x, X} at time t is given by 



Vix,X,t) = — — Q{x,X,t)P{x,X,t) 

j dXfP{xf,Xf,tf\xo, Xo,Q) 



where 



P{x,X,t) = P{x,X,t\xo,Xo,0) (27) 
Q{x,X,t) = j dXfP{xf,Xf,tf\x,X,t) (28) 

The coordinates Xf of the solvent are integrated over since the trajectories are not con- 
ditioned over the solvent molecules. 

Using the method described in the previous sections, the exact generalized Langevin 
equations satisfied by the coordinates are 



+ ril'\t) (29) 
+ vP{t) (30) 



dxi 


Di dU 


OXi 


dt 


ksT dxi 


dXi 


D2 dU 


1 ^^^^ 

' """^ dx. 


dt 


ksTdXi 



where Di and D2 are resp. the diffusion constants for protein and water molecules and the 
Gaussian noises 'qf'^\t) satisfy the relation 
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<r^p)(t) r/;/'^^(t') >= 2D,^,6,,6{t - t') 



.(1.2)/-/ 



(31) 



Consider eq.(28). Let us show that the additional force term vanishes. Indeed, 
Q{x,X,t) = J dXfP{xf,Xf,tf\x,X,t) 



(32) 



because of space and time translation invariance. Due to the integration over Xf in eq.(32) 



we see that Q{x^ X, t) does not depend on X, and thus the new drift term in (30) is absent. 
Therefore the exact equations for the conditional probability in presence of solvent are 

aing 



doc 2 


Di dU 


dt 


UbT dxi 


dXi 


D2 dU 


dt 


ksT dXi 



+ 2Di- 



dxi 



(33) 
(34) 



Using the Trotter approximation (18), these equations become (using vector notations) 

Xf — X 1 



dx 

dt tf-t 47 
dX _ _ D2 dU{x,X) 
dt ~ knT dX 



^{tf-t)V,V{x,X)+r]W{t) 



(35) 
(36) 



where the noises are Gaussian, correlated according to eq.(31). 

We thus conclude that in presence of the solvent, the protein is evolved through a modified 
Langevin equation with the effective potential V{x,X), whereas the solvent molecules are 
evolved according to the standard Langevin equation in presence of the original potential 
U{x,X). 

Extension of this method to the case of the native state (see previous section) is imme- 
diate. 



VI. EXAMPLE: THE QUARTIC DOUBLE- WELL 

We now illustrate the method on the example of barrier crossing in Id (quartic potential). 

u{x) = l{x'-ir 
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Figure 2: Potential U{x) (in black) and potential V{x) (in red). 

This potential has two minima at x = ±1, separated by a barrier of height 1/4. Note 
that at low enough temperature, the potential V{x) has two minima at points close to ±1 
and one minimum at x = (from eq.(|2])). Note that V{x) is much steeper than U{x) and 
thus more confining, around its minima. 

The model can be solved exactly by solving numerically the Fokker-Planck equation or 
by diagonalizing the Hamiltonian. All the examples are performed at low temperature 
T = 0.05, where the barrier height is equal to 5 in units of ksT and the Kramers relaxation 
time, given by the inverse of the smallest non-zero eigenvalue of H, is equal to tk = 366.39. 

On fig|3| we present a long trajectory {tf = 1000) obtained by solving the Langevin 
eq.([T| for a particle starting at Xq = — 1 at time 0. The general pattern described in the 
introduction can be easily checked: the particle stays in the left well for a time of the order 
of 550, then jumps very rapidly into the right well, where it stays for a time of the order of 
200, then jumps back to the left well where it stays again a time equal to about 250. 
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1000 



Figure 3: Full Langevin trajectory during time tf = 1000 with 2 transitions between the minima 



The two crossings times are very short, and we display an enlargement of the first tran- 
sition in fig|4j 

As can be seen, the crossing time for this specific trajectory is approximately tc ~ 2.5, 
much smaller than the Kramers time. 

In figj5} we plot two examples of two trajectories conditioned to cross the barrier during 
a time tf = 5. The trajectory in black is obtained by solving the exact bridge eq.(14) by 



computing exactly (using a spectral decomposition) the matrix element of the evolution 



operator, while the trajectory in red is obtained by solving the approximate eq.(19) with 



the exact same sequence of noise ri{t). In the left figure, the 2 trajectories are barely 
distinguishable, whereas the agreement is not as spectacular on the right figure. 



Next we look at some observables, obtained by averaging over many trajectories. 
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Figure 4: Enlargement of the first transition region 




Figure 5: Two sets (a) and (b) of exact trajectories (in black) and approximate trajectories (in red) 



In figj6l we plot: in black the exact average x{t) (obtained by a full expansion over the 



eigenstates of H), in red the average x{t) over 2000 trajectories obtained by solving eq.(19) 



and in blue, the average x{t) obtained by reweighting the trajectories according to eq.(22) . 
Plot (a) is obtained for = 2, plot (b) for = 5 and plot (c) for = 10. 
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Figure 6: Average position as a function of time for (a) = 2, (h) tf = 5, {c) tf = 10. Black curve: 
exact. Red curve: approximate. Blue curve: reweighted 

As expected, we see that the discrepancy between the exact (black) and the approximate 
(red) average x{t) increases with tf. For times shorter than the transition time rc, the 
agreement is excellent, whereas for t/ = 10 > re, the agreement is not as good. However, 
we see that the reweighting procedure, although not perfect, improves drastically the quality 
of the average for large tf. 

One of the main defects which appears in the approximate theory is the following: In 
the exact theory, the transition between the 2 minima can take place at any time between 
and tf. By contrast, it seems that in the approximate theory, the transition is driven by the 
final state and takes place only in the end of the trajectory. This effect remains negligible 
as long as tf < tq but becomes important for tf >Tc. We illustrate this problem in figjTjfor 
tf = 10. On the left figure, the exact and approximate trajectories make their transition 
in the last part of the time, whereas in the right figure, the real trajectory crosses in the 
beginning while the approximate trajectory still crosses in the last part. 

However, as we are interested quantitatively only in the region where the particle crosses 
the barrier, one can make long runs of approximate trajectories: They will not be good ap- 
proximations of the real trajectories, except in the end of the trajectory where the transition 
to the final state occurs. 

VII. CONCLUSION 

We have presented in this paper a novel method to generate paths following the Langevin 
overdamped dynamics, starting from an initial configuration and conditioned to end in a 
given final configuration (point or region of configuration space). We propose an approxima- 



17 



(b) 



(a) 



Figure 7: Two sets (a) and (b) of trajectories conditioned to cross the barrier. In black, exact 
trajectories and in red, approximate trajectories. 



tion which is vahd for small times. We have not been able to quantify how small should the 
time be, but the approximate dynamics seems to correctly reproduce the transition through 
a barrier. The approximate dynamics seems to have a tendency to confine the system in 
its initial configuration, and to allow for the transition only in the final stages. But this is 
not really a drawback since if we evolve approximately the system over long times, it will 
remain close to its initial condition, thus generating unreliable trajectories. However in the 
latest stages, the system will make its transition to the final state during a short time for 
which our approximation is reliable. 

One of the great advantages of this method is that all generated trajectories are statisti- 
cally independent. It is thus very easy to generate many of these trajectories using parallel 
computers. In addition, the trajectories can be reweighted to provide a faithful sample of 
the exact stochastic dynamics. Finally, this reweighting technique allows for the calcula- 
tion of the matrix element of the evolution operator, and thus allows for the generation of 
adequately sampled paths. 

The paths generated by our method can also be used either as initial paths to perform 
Monte Carlo transition path sampling, or as initial conditions for path minimization to 
determine Dominant Folding Paths. 

The method is as simple to implement as ordinary Langevin dynamics, and its application 
to simple models of protein folding is currently under way. 
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